Hexagonal arrangement of phospholipids in bilayer membranes
Chen Xiao-Wei, Yuan Ming-Xia, Guo Han, Zhu Zhi
Terahertz Technology Innovation Research Institute, Shanghai Key Laboratory of Modern Optical System, Terahertz Science Cooperative Innovation Center, School of Optical-Electrical Computer Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China

 

† Corresponding author. E-mail: zhuzhi@usst.edu.cn

Project supported by the National Natural Science Foundation of China (Grant No. 11904231), the National Key R&D Program of China (Grant Nos. 2018YFE0205501 and 2018YFB1801500), and Shanghai Sailing Program, China (Grant No. 19YF1434100).

Abstract

The phospholipid membrane plays a key role in myriad biological processes and phenomena, and the arrangement structure of membrane determines its function. However, the molecular arrangement structure of phospholipids in cell membranes is difficult to detect experimentally. On the basis of molecular dynamic simulations both in a non-destructive way and at native environment, we observed and confirmed that the phospholipids self-assemble to a hexagonal arrangement structure under physiological conditions. The underlying mechanism was revealed to be that there are hexagonal arrangement regions with a lower free energy around each lipid molecule. The findings potentially advance the understanding of biological functions of phospholipid bilayers.

1. Introduction

Phospholipid bilayers are the major building blocks of cell membranes in most living organisms, viruses,[1,2] and vesicles.[3,4] This bilayer has a fundamental bearing on various biological phenomena and processes, including the permeation of molecules,[5,6] drug delivery,[7,8] fertilization,[9,10] phagocytosis,[11,12] and the regulation of salt concentrations and pH values.[1315] In particular, the phospholipid bilayer participates directly or indirectly in the signaling process,[3,1618] including the transmission of nerve impulses[16] and immunological recognition,[17] as well as plays an important role in signaling or anchoring other molecules in cell membranes.[18]

The arrangement of structure determining the function is a key idea in biology.[19,20] Progresses have been made that the morphology of nanostructure determines the wettability, permeability, and adhesion nature of the nanostructure.[2123] For example, in 2018, we proposed a wettability change caused by the curvature vary of carbon-based and platinum-based surfaces.[24] Therefore, it is necessary to understand the arrangement information of the phospholipid in bilayer before its biological function, as well as the mechanism underlying the arrangement. Unfortunately, lipid bilayers are so thin and fragile that their arrangement structure is difficult to study with a traditional light microscope.[25] Fortunately, an increasing number of advanced technological means, i.e., x-ray crystallography (XRD), nuclear magnetic resonance (NMR), electronic microscopy (EM), fluorescence microscopy (FM), and atomic force microscopy (AFM), have been developed to detect the structure of the material in a destructive way or at a non-native environment. For example, experiments from EM[26,27] and AFM[28] suggest that the polar heads and the hydrophobic chains of lipids in the bilayer form hexagonal microcrystals. There is an urgent need to confirm the arrangement of lipid inside bilayer both in a non-destructive way and at native environment. Grazing-incidence scattering techniques[29] and theoretical calculations should be a promising solution. Besides, the mechanism underlying the arrangement of lipid in bilayer is call for elucidation.

Fortunately, classical molecular dynamic (MD) simulation is an effective mean to disclose the physics behind plentiful interesting phenomena on the basis of a micro perspective,[3037] and has also been widely used to study the characteristics of bilayers under native conditions.[3842] In 1982, Berendsen et al. as the pioneers applied the method of MD simulation to the representation of a realistic lipid bilayer, and found that this method can perfectly reproduce the experimental order parameters of the bilayer.[38] In 2002, Marrink et al. discovered the mechanical and electrical stress dependence of the pore formation and membrane rupture in phospholipid bilayers through MD simulations.[39] In 2013, Zhou et al. used MD simulations to reveal the underlying mechanism of the destruction of phospholipid membranes due to the toxicity of nanoparticles and nanosheets.[40] In 2015, Netz et al. quantitatively investigated the long-debated hydration repulsion of the biological membrane using MD simulations, and reviewed the physical mechanisms of the interaction between lipid membranes in an aqueous environment.[41] Recently, by means of MD simulations, Cremer et al. disclosed the complex nature of interactions between phospholipid bilayers with cations.[42]

In this article, we investigate the arrangement structure of lipid molecules in the bilayer under a native environment by MD simulations.[24,43] Since the major lipid component of the animal cell membrane is phosphatidylcholine (PC),[44] a bilayer consisting of the most common PC molecules, dimyristoyl-phosphatidylcholine (DMPC),[45] was studied. During several-microsecond simulations, we observed and further confirmed a hexagonal arrangement structure of phospholipids in the cell membranes at the thermodynamic equilibrium state. The underlying mechanism was revealed that there are hexagonal arrangements of the regions with a lower free energy for the free energy surface of each lipid molecule. Our finding thus advances the understanding of lipid arrangement structure in cell membrane as well as the biological function of the bio-membrane system.

2. Method

A bilayer system composed of 512 DMPC molecules, 21574 water molecules, and 0.15 mol/L NaCl was studied. The DMPC molecules and ions were characterized by the CHARMM 27 force field.[46] Water molecules were modeled by SPC/E.[47] The monomolecular structure of DMPC is shown in Fig. 1(a). Since carbon-21 (C21) is the center of two long hydrophobic tails of the DMPC molecule, we characterized the atom as the center of the DMPC molecule for convenience. The phospholipids bilayer was placed in the center of a box with an initial size of LX = 10.58 nm, LY = 12.85 nm, and LZ = 9.12 nm (see Figs. 1(b) and 1(c)). All simulations were performed with NPT ensemble using the software GROMACS 5.1.2.[48] In the process of the self-assembly of the phospholipid bilayer, the main simulation was maintained at 300 K and 1 bar using a Nose–Hoover thermostat[49,50] and the pressure coupling method presented by Parrinello–Rahman,[51,52] respectively. Periodic boundary conditions were applied in all directions. The particle-mesh Ewald (PME)[53] method was used to treat the long-range electrostatic interactions. The MD simulations were performed for 2 μs with a time step of 1 fs, and all the trajectories during the 2 μs simulation were collected for analysis with an interval of 1 ps.

Fig. 1. Schematic architecture of the simulated phospholipid bilayer system. (a) Monomolecular structure of DMPC. The red ball represents the C21 atom of a DMPC molecule, which is applied for labeling the location of the entire molecule. (b) Top view of the bilayer consisting of 512 DMPC molecules, 256 molecules per layer. (c) Side view of the simulated system, including the DMPC-containing bilayer, Na+ (blue balls), Cl (yellow balls), and water (light blue).
3. Results and discussion

The simulation results show that the lipid bilayer system composed of 512 DMPC molecules reached a dynamic equilibrium state on a microsecond time scale. In order to remove the influence of the artificial bilayer structure, we performed a trapezoidal heating and annealing process in the forepart simulation for 700 ns. Specifically, the simulated system was linearly heated from a temperature (T) of 300 K to 330 K in the previous 200 ns, maintained at 330 K in the following 300 ns, and linearly annealed from 330 K to 300 K in the last 200 ns. After heating and annealing, the self-assembly process of the bilayer was investigated for 1.3 μs under ambient conditions, which we denoted as the main simulation. The T of the simulated system was in accord with our expectations of the heating and annealing processes, and it kept equilibrium at 300 K in the main simulation (see Fig. 2(a)). Accompanying the control of the system temperature, the total energy of the system (E) and the average area per DMPC molecule (ρ) were changed correspondingly. Figure 2(a) also shows that E increases and decreases as T increases and decreases, respectively. This change is attributed to the kinetic energy of the system, an important component of E, which has a linear relationship with T summarized as the equipartition theorem. During the self-assembly process in the main simulation, E decreased lightly and finally converged to a value of −5.8 × 105 kJ/mol. Furthermore, ρ converged and fluctuated about a value of 0.45 ± 0.1 nm−2 during the process. Both the E and ρ were convergent on the microsecond time scale, so we deduced that the simulated bilayer system reached a dynamic equilibrium state. Intriguingly, observing the trajectory of the bilayer at the dynamic equilibrium state, we found a common phenomenon of six DMPC molecules regularly rounding a DMPC molecule. Figure 2(b) shows a typical one-sided conformation of the top view of a bilayer, where a series of paralleled hexagons were applied to cover the neighboring molecules of the DMPC molecule. The majority of the hexagons cover six DMPC molecules, which implies that phospholipid molecules in each layer of the bilayer may self-assemble into a hexagonal arrangement structure.

Fig. 2. Self-assembly dynamics of the bilayer system and the typical conformation of the bilayer at a dynamic equilibrium state. (a) Black, blue, and orange curves represent the time dependence of the system temperature (T), system total energy (E), and average area (ρ) per DMPC, respectively. (b) Typical locations of DMPC molecules from a top view of the bilayer. Red spheres represent the location of DMPC molecules. The endpoints of the black hexagons express the hexagonal arrangements of molecules around a DMPC molecule, where the hexagons are parallel to each other.

We further studied the distribution of neighboring DMPC molecules around a DMPC molecule. A vector Vij pointed from the center of the i-th DMPC molecule to the j-th one (see Fig. 3(a)). By means of the vector, the density distribution of neighboring molecules σ(r) = 〈Ni(r)〉/S with respect to a distance r = |Vij| and the angle distribution of neighboring molecules p(θ) = 〈Ni(θ)〉/∫N(θ)dθ with respect to an angle θ = acos {VijX/(|Vij||X|)} − αij ⋅ 360° could be obtained. Herein, Ni(r) and Ni(θ) are the numbers of neighboring DMPC molecules around the i-th molecule with a distance r and an angle θ in the same layer, respectively. S is the area. X is the unit vector of the x-axis for the simulation box, and αij is the switch constant of θ with a value of 0 or 1 when the projection value of Vij on the y-axis is positive or negative, respectively. The angular bracket denotes the average over all DMPC molecules. Figure 3(b) shows that σ (r) has two peaks under both high (330 K) and room (300 K) temperatures, which means that each DMPC molecule has two shells consisting of neighboring DMPC molecules inside the monolayer. When | Vij | is less than 0.6 nm, the j-th molecule belongs to the nearest neighbor shell of the i-th molecule. The second neighbor shell of the i-th molecule labels the j-th molecule whose | Vij | is larger than 0.6 nm and smaller than 1.2 nm. The results of p(θ) are shown in Fig. 3(c). Under the room temperature, p(θ) has six peaks for the nearest neighbor shell of a DMPC molecule. The peaks are located at the angles of −148°, −85°, −26°, 32°, 95°, and 154°. Likewise, there are also six peaks of p(θ) for the second neighbor shell with θ at −178°, −115°, −56°, 4°, 65°, and 124°. Remarkably, the peaks of p(θ) have an almost equal interval of ∼ 60° for both the nearest and the second neighbor shells of a DMPC molecule. Moreover, the angle displacements between the peaks of p(θ) for the nearest shell and the peaks for the second neighbor shell are approximately 30°. These results confirmed that the lipid bilayer self-assembles into a hexagonal arrangement structure under room temperature. For the case of high temperature, p(θ) of the nearest neighbor shell has four peaks with θ at −160°, −20°, 20°, and 160°; the corresponding four peaks of the second neighbor shell appear at angles of −108°, −62°, 73°, and 118°. Therefore, the hexagonal arrangement structure of lipids in the cell membrane exists at room temperature rather than at high temperature.

Fig. 3. Distribution of neighboring lipid molecules around a DMPC molecule. (a) Schematic diagram of molecular locations under a hexagonal arrangement structure. Red balls represent the locations of molecules; Vij is the vector from molecule i to j. (b) Radial distribution density of neighboring molecules around a DMPC molecule, black and red lines show the cases at 330 K and 300 K, respectively. (c) Probability distribution (p) of neighboring molecules around a DMPC molecule with respect to the angle θ between Vij and the x-axis of the simulation box. The upper and lower panels represent the cases where Vij points from a DMPC molecule to the nearest neighbor and the second-nearest neighbor, respectively. Gray circles and red squares correspond to the conditions with system temperatures at 330 K and 300 K, respectively.

Free energy landscape allows the self-assembly process of bio-macromolecules to be described and visualized in a meaningful manner.[54] We thus investigated the relative free energy landscape of a DMPC molecule. Two parameters rX = |Vij| ⋅ cosθ and rY = |Vij| ⋅ sinθ were applied to explore the free energy at the position (rX, rY) around a DMPC molecule. The free energy fluctuation can be estimated using the following equation:[55,56]

where P(rX, rY) is the probability of neighboring DMPC molecules falling into a grid which has an area of 0.5 Å × 0.5 Å and centers around the point (rX, rY); NA is the Avogadro constant; kB is the Boltzmann constant; and T is the average temperature of the simulated system. Figure 4(a) shows that there are six regions with low free energy around each DMPC molecule. The regions are in line with the nearest neighbor shells of any one molecule in the structure with a hexagonal arrangement consisting of lipids. Likewise, another six regions with a lower free energy also exist farther away from the molecule, corresponding to the second neighbor shells of the hexagonal arrangement structure. Each DMPC molecule has a similar free energy landscape, so phospholipid molecules in the bilayer are likely to self-assemble concertedly into a hexagonal arrangement structure as 3-dimensionally shown in Fig. 4(b).

Fig. 4. Free energy landscape of lipid molecule in bilayer. (a) Landscape of the free energy fluctuation (F) around a DMPC molecule with respect to parameters rX and rY. The red and blue colors correspond to the free energy fluctuation with lower and higher values, respectively. The unit of free energy fluctuation is kJ/mol. The solid and dashed hexagons represent the nearest and second neighbor regions of a DMPC molecule, respectively. (b) 3-dimensional arrangement sketch of lipids in bilayer membrane, where a cylindrical lipid molecule consists of a red ball and two gray curves.
4. Conclusion and perspectives

To summarize, MD simulation was employed to investigate the arrangement structure information of phospholipids in cell membranes under native conditions. On the basis of the simulation, we found that lipids in cell membranes can self-assemble into a hexagonal arrangement structure under physiological conditions. The underlying mechanism is attributed to the free energy surface of each DMPC that exists in several regions with hexagonal arrangements, where the regions possess a lower free energy.

It is worth noting that phosphatidylcholine is a strong polar molecule, which is tunable by an electromagnetic signal. Recent research disclosed that biological membranes might play a role in the transmission waveguide of light information in vivo.[5759] Advances in lipid arrangement structure in biological membranes may provide an understanding of the recording, duplication, and transmission of information in biology.

Reference
[1] Heimburg T Jackson A D 2005 Proc. Natl. Acad. Sci. USA 102 9790
[2] Edidin M 2003 Nat. Rev. Mol. Cell. Bio. 4 414
[3] Langton M J Scriven L M Williams N H Hunter C A 2017 J. Am. Chem. Soc. 139 15768
[4] Bhaskara R M Linker S M Vogele M Kofinger J Hummer G 2017 Acs Nano 11 1273
[5] Yang K Ma Y Q 2010 Nat. Nanotechnol. 5 579
[6] Khuntawee W Wolschann P Rungrotmongkol T Wong-ekkabut J Hannongbua S 2015 J. Chem. Inf. Model. 55 1894
[7] Prausnitz M R Mitragotri S Langer R 2004 Nat. Rev. Drug Discov. 3 115
[8] Prausnitz M R Langer R 2008 Nat. Biotechnol. 26 1261
[9] Suzuki Y Nagai K H Zinchenko A Hamada T 2017 Langmuir 33 2671
[10] Chen E H Olson E N 2005 Science 308 369
[11] Champion J A Mitragotri S 2006 Proc. Natl. Acad. Sci. USA 103 4930
[12] Arandjelovic S Ravichandran K S 2015 Nat. Immunol. 16 907
[13] Cheng Y L Bushby R J Evans S D Knowles P F Miles R E Ogier S D 2001 Langmuir 17 1240
[14] Cady S D Schmidt-Rohr K Wang J Soto C S DeGrado W F Hong M 2010 Nature 463 689
[15] Sharma M Yi M G Dong H Qin H J Peterson E Busath D D Zhou H X Cross T A 2010 Science 330 509
[16] Stevens C F 1993 Cell 72 55
[17] Simons K Toomre D 2000 Nat Rev Mol. Cell Bio. 1 31
[18] Simons K Ikonen E 1997 Nature 387 569
[19] Liu J Barry C E Besra G S Nikaido H 1996 J. Biol. Chem. 271 29545
[20] Beattie M E Veatch S L Stottrup B L Keller S L 2005 Biophys J. 89 1760
[21] Buten C Kortekaas L Ravoo B J 2019 Adv. Mater. 1904957
[22] Perez-Gonzalez C Alert R Blanch-Mercader C Gomez-Gonzalez M Kolodziej T Bazellieres E Casademunt J Trepat X 2019 Nat. Phys. 15 79
[23] Su B Tian Y Jiang L 2016 J. Am. Chem. Soc. 138 1727
[24] Zhu Z Guo H K Jiang X K Chen Y C Song B Zhu Y M Zhuang S L 2018 J. Phys. Chem. Lett. 9 2346
[25] Tieleman D P Marrink S J Berendsen H J 1997 Biochim. Biophys. Acta 1331 235
[26] Pilgram G S Vissers D C van der Meulen H Koerten H K Pavel S Lavrijsen S P Bouwstra J A 2001 J. Invest. Dermatol. 117 710
[27] Hui S Parsons D Cowden M 1974 Proc. Natl. Acad. Sci. USA 71 5068
[28] Muscatello U Alessandrini A Valdré G Vannini V Valdré U 2000 Biochem. Biophys. Res. Commun. 270 448
[29] Jaksch S Gutberlet T Müller-Buschbaum P 2019 Curr. Opin. Colloid & Interface Sci. 42 73
[30] Yu X M Qi C H Wang C L 2018 Chin. Phys. B 27 060101
[31] Qin J Y Geng Y Z Lu G Ji Q Fang H P 2018 Chin. Phys. B 27 028704
[32] Fang G Sheng N Jin T Xu Y S Sun H Yao J Zhuang W Fang H P 2018 Chin. Phys. B 27 030505
[33] Dou Q T Zuo G H Fang H P 2012 Chin. Phys. Lett. 29 068701
[34] Li Z C Duan L L Feng G Q Zhang Q G 2015 Chin. Phys. Lett. 32 118701
[35] Zhao L Tu Y S Wang C L Fang H P 2016 Chin. Phys. Lett. 33 038201
[36] Jiang Z L Chen P R Zhong W R Ai B Q Shao Z G 2018 Acta Phys. Sin. 67 226601 in Chinese
[37] Zhu Z Chang C Shu Y S Song B 2020 J. Phys. Chem. Lett. 11 256
[38] Van der Ploeg P Berendsen H J C 1982 J. Chem. Phys. 76 3271
[39] Tieleman D P Leontiadou H Mark A E Marrink S J 2003 J. Am. Chem. Soc. 125 6382
[40] Tu Y S Lv M Xiu P Huynh T Zhang M Castelli M Liu Z R Huang Q Fan C H Fang H P Zhou R H 2013 Nat. Nanotechnol. 8 594
[41] Schlaich A Kowalik B Kanduc M Schneck E Netz R R 2015 Physica A 418 105
[42] Bilkova E Pleskot R Rissanen S Sun S Czogalla A Cwiklik L Rog T Vattulainen I Cremer P S Jungwirth P Coskun U 2017 J. Am. Chem. Soc. 139 4019
[43] MacKerell A D Bashford D Bellott M Dunbrack R L Evanseck J D Field M J Fischer S Gao J Guo H Ha S Joseph-McCarthy D Kuchnir L Kuczera K Lau F T K Mattos C Michnick S Ngo T Nguyen D T Prodhom B Reiher W E Roux B Schlenkrich M Smith J C Stote R Straub J Watanabe M Wiorkiewicz-Kuczera J Yin D Karplus M 1998 J. Phys. Chem. B 102 3586
[44] Uran S Larsen A Jacobsen P B Skotland T 2001 J. Chromatogr B 758 265
[45] Ravula T Ramadugu S K Di Mauro G Ramamoorthy A 2017 Angew. Chem. Int. Edit. 56 11466
[46] Feller S E MacKerell A D 2000 J. Phys. Chem. B 104 7510
[47] Berendsen H Grigera J Straatsma T 1987 J. Phys. Chem. 91 6269
[48] Abraham M J Murtola T Schulz R Pall S Smith J C Hess B Lindahl E 2015 Software X 1 19
[49] Nose S 2002 Mol. Phys. 100 191
[50] Hoover W G 1985 Phys. Rev. A Gen. Phys. 31 1695
[51] Nose S Klein M L 1983 Mol. Phys. 50 1055
[52] Parrinello M Rahman A 1981 J. Appl. Phys. 52 7182
[53] Darden T York D Pedersen L 1993 J. Chem. Phys. 98 10089
[54] Dinner A R Sali A Smith L J Dobson C M Karplus M 2000 Trends Biochem. Sci. 25 331
[55] Hummer G Rasaiah J C Noworyta J P 2001 Nature 414 188
[56] Qi W P Song B Lei X L Wang C L Fang H P 2011 Biochemistry 50 9628
[57] Kumar S Boone K Tuszynski J Barclay P Simon C 2016 Sci. Rep. 6 36508
[58] Kwon J Kim M Park H Kang B M Jo Y Kim J H James O Yun S H Kim S G Suh M Choi M 2017 Nat. Commun. 8 1832
[59] Song B Shu Y S 2020 Nano Res. 13 38